function b_AFAD=Clogit_AFAD(Y_AFAD,X_AFAD,Z_AFAD)

choice=2:13;
ns=numel(choice);

Y=sortrows(Y_AFAD,[1,2]);
Y=Y(:,3:end);
Y=Y-1;  %!!!### UPDATED on 06/03/2020: Y MUST be from 1 to J!!!

X=sortrows(X_AFAD,[1,2]);
X=[ones(size(X,1),1),X(:,3:end)];

Z_AFAD=sortrows(Z_AFAD,[1,2,3]);
Z_AFAD(Z_AFAD(:,3)==1,:)=[];
Z=zeros(size(Y,1),size(Z_AFAD,2),ns);
for i=4:size(Z_AFAD,2)
    Z(:,i,:)=reshape(Z_AFAD(:,i),ns,[])';
end
Z(:,1:3,:)=[];


% X coefficients: constant bdtot perc_mcr perc_mcd hhi totpop65  
b(:,3) = [-3.490; 0.00118; 1.323; 1.389; 0.605; 0.241];
b(:,4) = [4.656; -0.00460; 3.264; 2.543; -1.273; -0.543];
b(:,5) = [2.636; -0.00512; 2.653; -2.312; -0.366; -0.360];
b(:,6) = [-3.855; 0.00220; 2.672; 1.073; -1.098; 0.157];
b(:,7) = [2.084; 0.00244; 1.535; -1.984; -1.973; -0.200];
b(:,8) = [-4.135; 0.00185; 0.735; -2.861; 0.863; 0.285];
b(:,9) = [4.030; -0.00554; 1.707; 1.581; -0.612; -0.420];
b(:,10) = [2.272; 0.000658; 2.895; 3.560; -2.360; -0.329];
b(:,11) = [-1.597; 0.000471; 0.273; 0.618; 0.698; 0.0651];
b(:,12) = [1.165; -0.0000457; 1.931; 1.543; -0.372; -0.151];
b(:,13) = [-2.909; 0.00110; 3.034; 2.352; 0.188; 0.0109];

% Z coefficients: vsysh vlag
bz=[1.916; 4.056];

bOri=[reshape(b(:,3:end),[],1);bz];
%bAns=zeros(size(bOri));


options=optimset('MaxFunEvals',200000000000,'MaxIter',1500000000000,'TolX',1e-16,'Tolfun',1e-16,'GradObj','on','DerivativeCheck','off','FinDiffType','central');
startval = bOri;
baseAlt=1;
restrMat=[];

b_AFAD = fminunc('clogit',startval,options,restrMat,Y,X,Z,baseAlt);








